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Abstract 

In this paper, we show how classical statistical field theory techniques 
can be used to efficiently perform the numerical evaluation of the non- 
perturbative Schwinger mechanism of particle production by quantum 
tunneling. In some approximation, we also consider the back-reaction of 
the produced particles on the external field, as well as the self-interactions 
of the produced particles. 

1 Introduction 

The Schwinger mechanism [1] (see [2] for a comprehensive review) is a phe- 
nomenon by which charged particles, e.g. electron-positron pairs, are pro- 
duced spontaneously from an external electrical field. This phenomenon is 
non-perturbative since pairs can be produced even from a static electrical field, 
something which is forbidden at any finite order of perturbation theory by sim- 
ple kincmatical arguments. It is also a purely quantum phenomenon, whose 
probability goes to zero in the classical limit h — > 0. Loosely speaking, e+e~ 
vacuum fluctuations are promoted to on-shell real particles by picking energy 
to the electrical field - which can be viewed as a kind of quantum tunneling 
process. 

In Quantum Electrodynamics, the probability for particle production by the 
Schwinger mechanism is of the order of exp(— irm 2 /eE), for particles of mass 
m and electrical charge e, in a field E. For fields one may realistically create 
in experiments, and taking even the lightest charged particle, the electron, this 
probability is so small (mostly due to the fact that the coupling constant e is 
small) that this phenomenon has remained elusive in all laboratory experiments 
so far (the typical electrical field necessary to make the production of an electron- 
positron pair likely is of the order of E ~ m 2 /e ~ 10 18 V/m). 



The subject of pair production by the Schwinger mechanism is also relevant 
in the context of Quantum Chromodynamics and strong interactions, since the 
strong coupling constant g is much larger. It is for instance an important in- 
gredient in hadronization models such as the Lund string model [3], where the 
breaking of a "string" made of a color electrical field into quark-antiquark pairs 
leads to meson production. It is also an ingredient in several phenomenological 
models of heavy ion collisions, e.g. [4-9]. 

It may also be a relevant mechanism of particle production in the Color Glass 
Condensate framework (see [10-14]), that is commonly used in the description of 
the first stages of hadronic or nuclear collisions at high energy. In this effective 
theory, the fast partons -mostly gluons at high energy- of the two colliding 
projectiles act as a static classical color source. The gluon occupation number, 
and therefore also this color source, increases with energy. Eventually, when the 
gluon occupation becomes of order of the inverse strong coupling 1/g 2 , nonlinear 
effects that tame this growth become important - an effect known as gluon 
saturation [15-17]. In this regime, the color source corresponding to the fast 
partons is of order 1/g, and therefore it creates fields that are themselves of order 
1/g. The probability of pair creation by such a strong field is not suppressed 
since gE can be large, unlike in QED. In [18], it has been argued that the CGC 
framework at next-to-leading order (1-loop) includes the contribution of the 
Schwinger mechanism to particle production. 

The Color Glass Condensate provides a semi-classical description of the un- 
derlying dynamics: at leading order (tree level), observables are computed by 
solving classical field equations of motion. This power counting is justified by 
the large occupation numbers and large fields that characterize the saturation 
regime, that allow one to neglect the non-commuting nature of the quantum 
fields and treat them as classical. The color fields obtained are leading order 
-the Glasma [19,14]- are non-pcrturbatively large, of order 1/g, and can thus 
lead to a large pair production. At next-to-leading order (1-loop), one can for- 
mulate observables in terms of small perturbations of this classical field [20,21], 
that obey linearized (and still classical) equations of motion. Equivalcntly, these 
1-loop corrections can be calculated, and resummed, by computing a classical 
path integral where one sums over a Gaussian ensemble of initial conditions for 
the classical field encountered at leading order. This approach, sometimes called 
classical statistical field theory, has been employed in a number of problems in 
cosmology [22-24], cold atom physics [25,26], and more recently in computations 
related to the thermalization of the quark-gluon plasma in heavy ion collisions 
[27-32]. 

The traditional method of computing the Schwinger mechanism has been 
to obtain it from the imaginary part of the Heisenberg-Euler Lagrangian [33], 
which in more modern quantum field theory language corresponds to calculating 
the imaginary part of the 1-loop effective action. It has also been derived in the 
framework of kinetic theory [34-36]. In Refs. [37-39], the Schwinger mechanism 
was computed from the Bogoliubov transformation that maps the creation- 
annihilation operators at t = +oo onto those at t = — oo. This requires that one 
solves the linearized equations of motion in order to obtain the time evolution 
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of the mode functions (i.e. modes that start as plane waves in the remote past, 
and are distorted by their propagation over the external field). 

Motivated by the applications of the classical statistical method to the CGC 
framework, we show in the present paper how the Schwinger mechanism can 
be calculated in classical statistical field theory, despite being an intrinsically 
quantum phenomenon. In order to keep the formalism as light as possible, we 
consider scalar electrodynamics instead of QCD. Note that a similar approach 
has been used in the case of fcrmion production in ref. [40-42] , following an idea 
of ref. [43] to simulate fermions efficiently on a lattice. 

In the section 2, we describe the model and discuss two methods of cal- 
culating the spectrum of produced particles at leading order; first in a rather 
standard quantum field theory formulation, and secondly as a classical path 
integral. In the section 3, we describe the lattice formulation of this calculation 
in order to evaluate the particle spectrum numerically, and then we compare 
the results of the classical statistical approach to the results obtained by the 
direct calculation of the I-loop diagram, in order to show that they are indeed 
equivalent. We improve this calculation in the section 4 in order to include the 
back-reaction of the produced particle pairs on the electrical field. Indeed, this 
screening effect is crucial for proper energy conservation. The self-interactions 
among the produced particles, that arc crucial for the eventual thermalization 
of the system, are considered in the section 5. Since the issue of thermalization 
in classical statistical field theory has been addressed elsewhere, we focus here 
on the issue of mass renormalization, which plays a crucial role in the Schwinger 
mechanism due to its extreme sensitivity on the mass of the particles being pro- 
duced. Finally, the section 6 contains some concluding remarks, while details 
about the mass renormalization are relegated to the appendix A. 

2 Single inclusive spectrum at leading order 
2.1 Scalar QED model 

Let us consider the case of a complex scalar field (f> with U(l) symmetry, min- 
imally coupled to an Abelian vector field A 1 *. The vector field may be coupled 
to an external source J M that drives it to a non-perturbatively large value. The 
Lagrangian of this model is 

C = -\F^F» V + (D^)(D^)* - m 2 <j>*<f> - *W) + 

= d^A v - d u A tl , = - ieA^ , (I) 

where e is the electrical charge of the scalar particles described by the field <f). 
We have not specified for now the self-interaction potential V of the scalar field, 
except for the fact that it depends only on the U(l) invariant 00*. A typical 
example of such a potential would be a quartic interaction, 

V(<M>*) = (2) 
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where A sets the strength of the self-interactions. 

The external source J^ xt can produce a non-trivial gauge potential, which 
in turn may produce scalar particles. Assuming that the initial state of the 
system is the vacuum, the inclusive spectrum 1 of scalar particles is given by the 
following formula in terms of the 2-point correlation function of the held <j>, 

dN x 1 



, ? , w>E.. I d4xd4y eiHX ~ V) (□* + ™ 2 )0=>v + ™ 2 ) <0 to |^(ar)^(y)|0 i »> , 

"" ' [' : (3) 
where E 2 = p 2 + m 2 . This expression is simply the Lehmann-Symanzik- 

Zimmcrmann reduction formula for the expectation value of the number op- 
erator a)(p)a(p). Although this form of the reduction formula involves 4- 
dimensional integrals over the entire space-time, it can also be written in terms 
of purely spatial integrals thanks to the identity 

J d A x (D x + m 2 )4>{x) = J d 3 x [<j>(x°, x) - iE p cf>(x°, x)] - ( 4 ) 

where L4(x )]„ = A(b) — A(a) and the dot denotes a time derivative. Moreover, 
the lower boundary x° — > — oo docs not contribute since we assume that the 
initial state is empty 2 . Thus, eq. (3) can also be written as 3 

^1 - i im 1 f d 3 xd 3 v e - ip < x - y ) 

x (0 in | (t, x) + iE p d?* (t, x))(4>(t, y) - iE p cj)(t, y)) |0 in ) . (5) 

On may remove the limit t — > oo in this formula, and interpret its result as the 
particle spectrum at the time t. However, one has to keep in mind that this 
interpretation cannot be completely rigorous: strictly speaking, the particles 
need to be free and on-shell in order for their number to be a well defined 
concept, which takes an infinite time. 

In addition, this formula for the spectrum assumes that the gauge potential 
vanishes when t — > +oo. However, even if the electrical and magnetic fields are 
made to vanish in this limit, the gauge potential itself could be a non-zero pure 
gauge = d^x- In this case, we need to perform in eq. (5) the replacement 

e -ip x _^ e -ip x e iex(x) ^ 

(In words, we must gauge transform the free plane waves that are used in the 
Fourier decomposition of the fields) . Note that this replacement is also required 

1 This observable should not be confused with the probability Pi of producing exactly one 
particle-antiparticle pair, that would be obtained from the matrix element (0 O ut|</>(a;)|0i n ). 
The average number of produced particles (i.e. the integral over p of dNi/d?p) is related to 
the probabilities P n by Ni = nP n . Ni is usually easier to calculate than the individual 
P n , thanks to simplifications related to the completeness of the set of possible final states. 
The contribution of the lower boundary is in fact (0in|4 n (p)a in (p)|0i n ) = 0. 

3 At this stage, the object <j> is still an operator, and one should keep in mind that the 
commutator [0, <f>] is not zero. 
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in order to have a gauge invariant spectrum. In the particular case where the 
gauge potential A is spatially homogeneous, this substitution can also be written 
as 

e -ipx _^ p -i(p+eA)-x ^ 

and we recognize now the well known difference between the kinetic and canon- 
ical momenta 

Pcano = Pkin + eA ■ (8) 

By extension, we will also perform this substitution in order to define the particle 
spectrum in the presence of a homogeneous electrical field (i.e. when the gauge 
potential is not a pure gauge) . Again, one must keep in mind that the concept of 
particle number in the presence of a non-trivial background field is not rigorously 
defined, since the particles we are trying to count are not free particles. 



2.2 Spectrum at Leading Order 

A typical graph contributing to the spectrum is shown in the figure 1. The 




Figure 1: Generic contribution to the inclusive spectrum of produced scalar 
particles. The wavy lines represent the photons, while the solid black lines are 
scalars. The crosses terminating the photon lines represent the source J^ t . 

spectrum can be organized as a triple series expansion, in powers of the electro- 
magnetic coupling e, of the self-coupling A, and of the external source J^. t . 

When the external source J£ xt is large, possibly of order J cxt ~ 0(e~ 1 ), 
one may expect non-perturbative effects such as the Schwinger mechanism to 
become important. These contributions are usually non analytic in e, and it is 
necessary to resum all the powers of the combination e J^ xt in order to compute 
them. In the rest of the paper, we call Leading Order the result of this resum- 
mation, in which we include all powers of e accompanied by a power of J^ t , but 
no further corrections in e 2 or in A. Therefore, the graphs that contribute to the 
spectrum at leading order are considerably simpler (see the figure 2), since they 
have only one scalar loop dressed by insertions of a photon directly connected 
to the external source J/t*. 
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Figure 2: Graphs contributing to the inclusive spectrum at leading order. 

At this level of approximation, one can treat the gauge field attached to the 
scalar loop is a classical field A^ 1 that obeys the classical Maxwell's equations, 

V = JZ« , (9) 

and the 2-point correlation function that enters in the spectrum is expressible 
as follows 



d 3 q 

(Vjy + m 2 )<p ±q (x) = , lim <p Tq (x) = e ±i «' x , (10) 



(0 in \(j)*(x)<f>{y)\0 in ) Tn = I lo , z (p- q (x)ip +q (y) 



where X> M = — ieA^ is the covariant derivative constructed with the classi- 
cal background field. Note that the mode functions ip+ q and (p- q are mutual 
complex conjugates since the background gauge field is real (a requirement of 
unitarity). 

At this point, the calculation of the scalar particle yield at leading order 
has been recasted into a purely classical calculation, where one needs to solve 
the classical equation of motion {V^V* 1 + m 2 )(p = for each of the scalar 
mode functions f± q . In the special case of a static electrical field, eqs. (10) are 
equivalent to the classic result of Schwinger. Note that eqs. (10) do not imply 
that the particle spectrum at leading order is a classical quantity. Indeed, it is 
well known in the case of a static electrical field that the particles are produced 
by a quantum tunneling phenomenon, whose probability goes to zero if K — > 0. 
Instead, eqs. (10) should be viewed as an example of the general property that 
one-loop quantities can be written as quadratic forms in terms of fields that 
obey linearized classical equations of motion. 

A crucial aspect of this formulation is that these mode functions are specified 
by retarded boundary conditions in the remote past, where they behave as free 
plane waves. Using the identity (4) and the fact that in the limit x° — > — oo, we 
have (p- q (x) = e lq ' x , we simply have 

/ d 4 x e lp - x (a x +m 2 )(p_ q (x) = lim / d 3 x e lv ' x [tp- q (x°, x)-iE p <p_ q (x° , x)} . 

(11) 
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Thus, as illustrated in the figure 3, the physical interpretation of cqs. (10) is 
that in order to obtain the spectrum of produced particles, one should start 
in the remote past with negative energy plane waves (that are equivalent, by 
crossing symmetry, to having a positive energy antiparticle in the final state), 
that subsequently evolve over the classical gauge field A 1 *, and are projected 
at the final time on a positive energy plane wave. The momentum q of the 
incoming plane wave can be interpreted as the momentum of the antiparticle 
that must be produced along with the observed particle of momentum p, and 
therefore should be integrated out in order to obtain the particle spectrum. 




Figure 3: Diagrammatic representation of eqs. (10). 

This formulation provides an explicit numerical method for computing the 
yield at leading order for a general background field (for which it is not possible 
to solve analytically the equation of motion for <p± q ): discretize space and 
solve numerically the classical equation of motion for each momentum q of the 
reciprocal lattice. This is however a computationally expensive method, since 
this computation scales like the square of the number of lattice points. More 
precisely, the computation time scales as 

Nt x JV£ tt , (12) 

where Matt is the number of lattice points and Nt the number of time-steps 
used in solving the equations of motion. 

2.3 Reformulation as a Gaussian functional integral 

It is however possible to formulate the particle spectrum in an alternative way, 
that allows a more efficient computation based on a Monte-Carlo sampling in a 
functional space that we shall specify shortly. Since the evolution of the mode 
functions is causal, knowing them, as well as their first order time derivative 4 , 
on some Cauchy surface S (for instance, any surface of constant time x ) is 
sufficient to fully determine their subsequent evolution. To illustrate this, let us 
consider the following functional 

G xy [<po,ir } ee <p*(x)ip(y) 

(V^ + m 2 )<p = , <p(ta,x)=<po(x), (p(t ,x) = iro(x). (13) 

4 This is necessary because the equation of motion for tp contains second order time deriva- 
tives. 
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Gxyifa, To] is the product at the points x and y of the classical solution (p whose 
initial conditions at the time to are given by the functions ipo(x) , ir (x) . Let 
us now average this functional over a Gaussian distribution W[<£o,7!"o] 01 these 
initial conditions, specified by the following 2-point correlation functions : 

Wo(x)My)) = J Jfr^jf ( P-q(to,x)ip+ q (t ,y) 
f d 3 q 

(TT*(x)TT (y)) = J 3^ <p- q {t ,x)ip +q {t ,y) 
(M^)Mv)) = 0. (14) 
Then, it is easy to check that 



/ 



[DtpoDiro] W[(p Q ,ir } G xy [ip , ir ] 
lE q 

= ^<Oi„|0*(x)0(y) + 0(y)0*(x)|O in ) Lo . (15) 

In other words, this procedure gives almost the expectation value we need in 
order to compute the particle spectrum, except for the ordering between the two 
field operators. In the form (5) of the reduction formula, the two field operators 
are evaluated at equal times. Therefore, <j)*{x) and <f>{y) commute, but not <j>*(x) 
and <p(y). In fact, if we used the expectation value (15) in the LSZ formula, 
we would get the expectation value of the operator (a i (p)a(p) + a(p)a i (p))/2 
instead of a)(p)a(p). This discrepancy is easy to fix by using the canonical 
commutation relation between creation and annihilation operators, 

Hp), a\q)] = (2tt) 3 2E p S(p - q) . (16) 

One can obtain the leading order spectrum by calculating the expectation value 
(15), and then by subtracting V/2 particles where V is the volume 5 of the 
system. More precisely, 



d 3 p 



V 

'~2 



J [DipoDno] W[<po,7r ] F p [<p , tt ] , (17) 



where 



1 r ,. 2 



F p [^o} = {2n)32Ep \* P [<Po,*o 

®p[(Po,7To]= Km / d s x e w ' x [(f{x°, x) - iE p (p(x°, x)] 

(£> M X> M + m 2 ) ip = , ip(t ,x) = <po(x) , <p(t ,x) = tt {x) . (18) 



5 The volume arises from a 5(0) in momentum space. 
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It is easy to check that if we carry out this procedure in the vacuum (i.e. with 
2-V = m the equation of motion of the scalar field), we obtain dN\/d?p = 0. 
The subtraction of the term V/2 is crucial for that. 

This reformulation of the spectrum at leading order can be illustrated dia- 
grammatically as follows 6 , 



(19) 



W[tp ,*a] F p [tpa,ir ] 

The new representation amounts to slicing the scalar loop at a certain time to, 
and to assign to the functional F p [<po, tto] the time evolution for x° > to, and to 
the Gaussian distribution W{tp ,n ] the part of the evolution at x° < t . The 
Gaussian average over the fields fo,^o is the "glue" that reconstructs the loop, 
since it amounts to connecting pairwise the two open endpoints of F p [cpo, 7To]. 
Note that the choice of the time to is not important, since the left hand side does 
not depend on this choice. In fact, instead of slicing the loop at a fixed time 
t , the separation in two factors could have been made on any locally space-like 
surface. In practical implementations, it is best to perform the separation at a 
time to such that the distribution VF^Oi To] is easy to compute. 

Note that random elements of the Gaussian ensemble defined by eq. (14) 
can be generated by writing 

V°( X ) = / (1 \20T? [ C q<P+q{t0,x) +d q tp- q {to 1 x)} 



(2Tr) 2 2E q 
f d 3 q 

= J (2tt) 2 2E [° g V+vitO' x ) + d i y-q(*o, a)] , (20) 

where c q and d q are random complex Gaussian distributed numbers, whose 
2-point correlations are 

(c q c ql ) = (d q d* q ,) = (2nfE q 6(q - q') 

(c q c q ,) = (d q d q/ ) = {c q d q i) = {c q d*q,) = . (21) 



2.4 Numerical cost 

Eq. (17) may at first appear to be a drawback compared to eq. (10) since we 
have replaced an ordinary integral by a functional integration. However, let us 
assume that the mode functions ip± q are known (or at least easily computable) 

6 Although in this illustration we have put all the interactions with the background field in 
the functional F p [ipo, tto], this is not mandatory. If the time to is chosen larger than the time 
at which the background field is switched on, then the Gaussian distribution W^i^Oii'o] a lso 
depends on the background field. 
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at the initial time to. This is for instance the case when the background field 
A^ is vanishing for x° < to- Then, one can estimate the functional integral by 
doing a Monte-Carlo sampling of the Gaussian ensemble defined by eqs. (14), 
i.e. by generating random functions of the form (20) and by solving the equation 
of motion for each of these samples. The computational time of this method 
would be of the order of 

^samples X N latt X (N Utt + N t ) , (22) 

where iV samp i es is the number of samples used in the Monte-Carlo evaluation of 
the functional integral. Note that the first term, proportional to iV^ tt , is due 
to the fact that in eq. (20) there is a sum over q at each position x. This has 
to be repeated for each sample, but needs to be done only at the initial time 
io- This Monte-Carlo method of evaluating the particle spectrum is less costly 
than the direct method provided that 

^samples < Matt , ^samples < N t ■ (23) 

The first condition is obvious: if one uses a number of samples which is larger 
than the number of independent mode functions, then one would be better off 
using the direct method (since it would give the exact leading order answer, for 
a lesser computational effort). The second condition implies that the computa- 
tion is dominated by the resolution of the equations of motion rather than the 
evaluation of the initial conditions. 



3 Lattice numerical evaluation 

3.1 Lattice setup 

In the two formulations, (10) or (17), one needs to solve the linearized equation 
of motion iV^V^ + in 2 )ip = for the propagation of a scalar field on top of 
some background electromagnetic potential A^. There are only a few known 
examples of background fields for which this equation of motion can be solved 
analytically. For a generic background field, one can only solve this equation 
numerically. 

Actual computations are done by discretizing space on a lattice. We consider 
a finite box of volume L T X Ly X L z and divide it into a iV latt = N x x N y x N z 
lattice. Space points are labeled as 

x = (n x a x ,n y a y ,n z a x ) (n x = !,■■■ ,N X ; n y = 1, • • • , N y ; n z = !,-■■ ,N Z ) , 

(24) 

with lattice spacing a x = L x /N x , etc. We impose the periodic boundary condi- 
tions, e.g. (p(x, y, z) — ip(x + L Xl y, z), and the momenta are given by 

pi= 2wK (i = Xj y jZ ^ (25) 
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,0," 



(i = x,y,z) 



(26) 



where the ki are integers, taken in the range 

2 

where we have assumed that the lattice sizes Ni are even. 

On the lattice, differentiation with respect to space is replaced by finite 
differences. Let us introduce the forward difference 



Vf<p(x) = — [ip(x + aihi) - ip(x)} 



(27) 



where the vector hi is the displacement by one lattice spacing in the spatial 
direction i. Similarly, the backward difference is 



V 4 <p(x) = — [<p(x) - <p(x - aihi)] 



(28) 



By "integration by parts" , the forward difference is transformed to the backward 
difference : 



J2 Mf{x)]g{x) = -'£f{x)\VTg{x) 



(29) 



(Notice that there is no boundary term because of the periodic boundary con- 
dition.) In words, this means that V~ and Vf are mutually adjoint. Since it is 
desirable to have a self-adjoint Laplacian operator, it is convenient to define it 
by a mix of forward and backward derivatives 



(30) 



The discrete plane waves exp(ip ■ x) are eigenfunctions of this operator, with 
eigenvalues 



sin 2 (^) 



a? 



(31) 

aj 

t—x,y.z 

When considering scalar QED, the local U(l) gauge invariance can be pre 
served on the lattice by defining the forward covariant derivative as 

1 



D+<p(x) 

Under a gauge transformation 7 

ip{t,x) - 
A (t,x) - 
Ai(t,x) - 



> A ^ x) Lp{x + aihi) - ip(x) 

e le0 ^ip(t,x) 
A (t,x)-6{t,x) 
Mt,x)-\7+9(t,x) , 



(32) 



(33) 



7 There is some arbitrariness in how we discretize the gauge transformation law for Aj, 
since we could have chosen a backward derivative instead of the forward derivative. If we 
adopt this alternative choice, the forward covariant derivative must be changed into 

( ,ieo i A 4 (a.+o 4 A i ) v ,( a . + a . ft .) _ ^ 



D+ip(x) 



-[< 

a; L 
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Dfip(t, x) is transformed in the same way as ip(t, x) : 

D+<p(t,x) — ► e iee ^D+v(t,x) . (34) 

One can write a gauge invariant discretized Lagrangian density for the complex 
scalar fields as follows 

Anattcr = (A)0)* (A>0 " E ( D t4>Y ( D t<t>) ~ ™ 2 4* <f> - V (0*0) . (35) 

i— x,y ,z 

In deriving the discretized classical equation of motion, one should note that the 
forward covariant derivative is the adjoint of the backward covariant derivative 

D-ip(x) = — L(x) - e - ieaiA ^ x - a ^ip(x - aihi)] . (36) 

One obtains 

(Dl- J2 Dr D + +m 2 )^ + V'{^v)v = Q- (37) 

i— x,y,z 

(Here the equation is written with the self-interaction term, but in the evaluation 
of the spectrum at leading order we need only the linear part of the equation.) 
It is convenient to choose the temporal gauge A — 0, so that D = do in the 
equation of motion. 



3.2 Numerical results 

In order to demonstrate that the classical statistical simulation (CSS) can indeed 
describe the Schwinger mechanism at leading order, we consider a simple situ- 
ation which can be handled easily by the direct quantum field theory method. 
The self-coupling A is set to zero in this leading order computation, and we take 
a spatially homogeneous background electrical field in the z direction (and no 
magnetic field), that we switch on at the time x° = 0. In this case, thanks to 
the translational invariance of the problem, the direct evaluation of eq. (10) is 
not very expensive and will be used as a benchmark against which we compare 
the CSS results. 

The parameters of the lattice simulation are N x = N y — 32, N z — 256 (we 
need more lattice spacings in the direction of the electrical field) , corresponding 
to physical sizes L x = L y = 50, and L z = 30 respectively. The mass is taken to 
be m = 0.1, the electrical charge is e = 1 and the electrical field is switched on 8 
at x a = and its value is E = 1 (alternatively, one may view eE as arbitrary 
and all the other dimensionful quantities as being quoted in units of y/ eE, that 
has the dimension of a mass). 1024 field configurations were used in order to 
sample the Gaussian ensemble defined in eqs. (20). 

8 In the temporal gauge A = 0, this can be realized by the following gauge potential : 
A 1 = A 2 = , -d°A 3 = Ee(x°) , 
which can be realized by the external current J^ t = — 8^ [i ES(x°). 
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In the figure 4 we show the longitudinal momentum distribution 9 of the pro- 
duced scalar particles, at several times shortly after the electrical field has been 
switched on. The dots represent the result of the classical statistical approach 
and the solid lines are the direct QFT calculation. The agreement between the 
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Figure 4: Comparison of the p z spectrum between the classical statistical simu- 
lation and the direct 1-loop QFT calculation in the case of a constant electrical 
field. 

two approaches is very good, and the differences compatible with the expected 
statistical error given the number of samples used in the CSS approach. In 
particular, the intricate oscillatory pattern of this spectrum, which results from 
quantum interference phenomena, is well reproduced in the CSS method. This 
shows very concretely how 1-loop quantum effects can be reformulated in terms 
of purely classical objects. 

The time evolution of the p z spectrum is rather transparent: particles are 
produced with a small p z by quantum tunneling, and later they are acceler- 
ated in the direction of the electrical field; hence the expansion of the spec- 
trum towards larger (positive, because we are considering only particles -not 
antiparticles-) values of p z . One can indeed see on the figure 4 that this ex- 
pansion is linear in time, in good agreement with a constant acceleration eE 
(which is equal to 1 with our choice of units) of the particles in the +z direction. 
Similarly, the comparison of the p± spectra obtained in the two approaches, in 
the figure 5, show a good agreement within statistical errors. The shape of the 
transverse momentum spectrum is very different from that of the longitudinal 
momentum distribution. Roughly speaking, the produced particles originate 
from virtual pairs (i.e. vacuum fluctuations) that can have a momentum in any 
direction. Then, their longitudinal momentum p z increases linearly in time due 
to the electrical field, while their transverse momentum is not affected. 

9 The occupation number /(p) differs from the particle spectrum defined in eq. (17) by a 
factor of the volume, dNi/dPp = V /(p)/(2vr) 3 . 
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Figure 5: Comparison of the p± spectrum between the classical statistical simu- 
lation and the direct 1-loop QFT calculation in the case of a constant electrical 
field. 

To close this section, we also show the energy density and electrical current 
density carried by the produced particles, as a function of time, in the figure 6. 
The energy density of the produced particles increases steadily with time, and 




Figure 6: Energy density and electrical current density carried by the produced 
particles. In the left plot, the light blue band indicates the energy density of 
the electrical field. 

at some point it will overcome the energy density of the background electrical 
field. When this occurs, the simple calculation done in this section, where the 
back-reaction of the produced particles on the gauge field is neglected, becomes 
certainly insufficient. In particular, the total energy present in the system is 
not conserved at this level of approximation, because of the uninterrupted pro- 
duction of scalar particles, shown in the left plot of the figure 6. One could in 
fact use this as a criterion for deciding when the one loop approximation ceases 
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to be valid: it should be improved when the energy carried by the produced 
particles becomes of the same order of magnitude as the energy density carried 
by the electrical field. From the left plot of 6, we expect this breakdown to occur 
as early as t ~ 10. This is also corroborated by the behavior of the electrical 
current carried by the produced particles. Because this current acts as a source 
in the Maxwell's equations that control the gauge potential, it can alter the 
background electrical field when it becomes too large. 



4 Back-reaction effects 

In the previous section, we have seen that the plain 1-loop calculation of the 
spectrum of produced particles is bound to break down after some time, because 
this approximation violates energy conservation. What is missing is the back- 
reaction of the produced scalar particles on the electrical field, which has the 
effect of screening this field, which eventually will put an end to the production 
of particles [44,45]. 

Taking into account the back-reaction means that the source of the elec- 
tromagnetic field is not just the external source J[? xt , but also includes the 
contribution of the scalar field </> to the electromagnetic current, 

j (x) = ie (j>* (x)<j)(x) — ((>* (x)(f>(x) 

f(x) = ie[<p(x)D+<t>(x)-(D+<t>(x)y<t>(xj\ , (38) 

written here with lattice discrete spatial derivatives. In the gauge A° = 0, the 
Maxwell's equation for the vector potential read 

v*_A l = J e ° xt +i° 

A 1 + (v;vi - 5* j 'v+ • v_) a 3 = .r cxt + f . (39) 

The first of the two equations (39) is Gauss's law. It is automatically satisfied 
if A 1 obeys the second equation, and is conserved 10 , 

j° = V!_J\ (40) 

which, in turn, is the case if the scalar field (j) obeys the equation of mo- 
tion. Diagrammatically, solving simultaneously the equation of motion of the 
scalar field (j> and Maxwell's equations, starting from some initial condition 
<Po(x), tto(x) at the time to, amounts to resumming graphs such as the one rep- 
resented in the left part of the figure 7. All the graphs are made of a principal 
scalar line, to which the measured particle of momentum p is attached. To this 
line are attached a number of photons. These photons can either be attached to 
the external current J^t, or to the induced current, i.e. to another scalar line. 
These secondary scalar lines can themselves be decorated by photons, etc... 



10 This follows from the identity V ! 1(V + • V_) = (V_ • V+)Vi 
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Figure 7: Left: topologies included when solving simultaneously the equations 
of motion of the scalar field (j> and the Maxwell's equations with the induced 
current, from some given boundary conditions ifo^Q at the initial time to- 
Right: QFT topologies when we perform a Gaussian average of the left graph 
over the Gaussian ensemble (14), with j M replaced by the ensemble average 
in Maxwell's equations. 

Averaging the functional of (fo, 7i"o obtained by this procedure over the Gaus- 
sian distribution of initial conditions defined by (14) amounts to reconnecting 
pairwise all the hanging scalar lines in the graph of the figure 7. Although 
all the topologies obtained when doing this were indeed present in the original 
quantum field theoretical formulation of the particle spectrum, some of them are 
miscalculated by the classical statistical simulation. This because in the CSS 
all propagators are retarded, while they are Feynman propagators in the field 
theoretical calculation 11 . However, when we reconnect together the two scalars 
that enter in the same instance of the induced current, the classical statisti- 
cal approach reproduces exactly the QFT value of that loop. This amounts to 
replacing the induced current j 71 in Maxwell's equation by its ensemble average, 

r (f) ■ (41) 

By doing this, one obtains all the graphs such as the one represented in the right 
part of the figure 7, where the scalar loops represented in orange originate from 
the use of (j^) in the right hand side of Maxwell's equations. Although this is 
not manifest on this example, these scalar loops can themselves be dressed by 
an arbitrary number of photon insertions, whose other endpoint can either be 
the external current J^ t or another instance of the ensemble averaged induced 
source. 

Let us now show some numerical results that illustrate the effect of the back 
reaction. First, in the figure 8, we display the electrical current density (left) 
and the resulting electrical field (right). There in an inflection in the growth of 
the current at a time t ss 10, which roughly corresponds to the moment when 
the energy carried by the produced particles becomes comparable to the energy 
stored in the electrical field. Simultaneously, the electrical field decreases and 
even changes sign periodically, while the amplitude of its oscillations decrease 
to zero. Consequently, the production of particles slows down and effectively 

lx In other words, the CSS and the original quantum field theory differ by some commutators, 
as one may expect. 
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Figure 8: Time evolution of the electrical current density and of the resulting 
electrical field when the back-reaction is taken into account. 



stops after some time (when the probability of particle creation, of the order 

of exp(— irm 2 /eE), becomes too small). It is rather easy to understand why 
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Figure 9: Illustration of the polarization phenomenon that leads to the reduction 
of the electrical field when the back-reaction is taken into account. 

the induced electrical field is opposite in direction to the externally field, as 
illustrated in the figure 9. Indeed, after having been produced, the positive 
charges are accelerated in the direction of the external field and the negative 
charges in the opposite direction. This charge separation acts like capacitor 
plates, that create an induced field oriented from the positive to the negative 
charges. The induced electrical field thus counteracts the effect of the external 
field. 

In the left plot of the figure 10, we see that the expansion of the p z spec- 
trum towards the right is no longer linear in time, a direct consequence of the 
decreasing electrical field which is no longer providing a constant acceleration 
to the produced particles. In fact the plot on the right of the figure 10 shows 
that for slightly larger times, the shift towards the right of the p z spectrum 
comes to a halt, and is replaced by a shift to the left. Obviously, this happens 
when the electrical field has changed sign, around t « 30. After the p z spec- 
trum moves into the negative momentum region, its value undergoes a rapid 
increase. This is because particles which have been created earlier pass through 
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the zero-momentum region and stimulate the subsequent particle production 
(Bose-enhancement). At the same time, the p z spectrum starts to show an 
oscillatory pattern. This is due to an interference between the fields of the 
previously produced particles and of the newly produced particles [38]. It is 
remarkable that the CSS method can describe this intricate pattern of peaks, 
that are purely quantum effects. At even larger times, when the electrical field 
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Figure 10: Left: evolution at short times of the p z spectrum of the produced 
scalar particles. Right: evolution of the p z spectrum at intermediate times. 

has become small, the p z distribution becomes roughly centered around p z = 0, 
as one can see from the plot on the left of the figure 11. Because the elec- 
tric field has changed its sign for several times, the spectrum has gone through 
several stages of Bose-enhancement and is now much larger than its value at 
early times. Finally, one can also check that the inclusion of the back-reaction 
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Figure 11: Left: evolution of the p z spectrum over longer time scales. Right: 
decomposition of the energy of the system between fields and particles. 

effects cure the main problem of the plain 1-loop calculation, i.e. the energy 
conservation. In the plot on the right of the figure 11, we have represented the 
energy density carried by the produced particles, the energy density carried by 
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the electromagnetic fields, and the sum of the two contributions. We observe a 
transfer into particles of the energy initially stored in the electrical field, while 
the total energy remains constant. 

5 Self-interactions and mass renormalization 

So far, all the results we have shown have neglected the self-interactions of 
the scalar fields. This means that after the external electrical field has been 
"neutralized" by the produced particles, their distribution is frozen and ceases 
to evolve. In particular, it has no way to thermalize. The effect of these self- 
interactions has been studied in the classical statistical framework employed 
here, where it is rather straightforward to take into account since it just amounts 
to adding a non-linear term in the equation of motion of the classical scalar fields, 

(Dl- ]T D^D* +m 2 )<p+±{<p*<p)<p = Q, (42) 

written here for a quartic self-interaction. In the slightly simpler example of a 
real scalar field theory, it has been shown that these nonlinearities lead to the 
isotropization and thermalization 12 of the momentum distribution of the parti- 
cles. Diagrammatically, including the self-interactions in the classical equation 




Figure 12: Left: topologies included when solving simultaneously the equations 
of motion of the scalar field <j> with self-interactions and the Maxwell's equations 
with the induced current, from some given boundary conditions (fo,TTo at the 
initial time to. Right: topologies obtained after the Gaussian average over the 
initial conditions ipo, ttq, with (j 71 ) as the source term in Maxwell's equations. 

of motion of the scalar fields changes the figure 7 into the figure 12. In particu- 
lar, the Gaussian average over the initial conditions for the scalar field can now 
produce loop corrections such as (but not only) tadpoles. 

Our goal in this section in not to reproduce previous results on thermalization 
in classical statistical simulations, but to stress the complication due to mass 
renormalization 13 , which is crucial when dealing with a tunneling phenomenon 

12 Since it is a semi-classical approximation, the asymptotic spectrum obtained in the classi- 
cal statistical framework is not a full fledged Bose-Einstcin distribution but only its soft part 

/Op) = t/e p . 

13 One can find other discussions of renormalization in classical approaches in refs. [46,47]. 
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such as the Schwinger mechanism. The main issue is that the probability of 
particle production by quantum tunneling is extremely sensitive to the value 
of the mass of the scalar particles, since its square enters in the exponential 
cxp(— irm 2 /eE). However, when we include the nonlinear term in the equation 
of motion of the scalar field and we average over its initial conditions, we resum 
some loop corrections that are ultraviolet divergent. In a lattice simulation, they 
are regularized by the lattice spacing, but provide a potentially large renormal- 
ization of the mass that will alter significantly the production of particles by 
the Schwinger mechanism. The worst offenders are the tadpole corrections, that 
have a quadratic dependence on the inverse lattice spacing. 
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Figure 13: p z spectrum of the produced scalar particles at t = 10, in the non 
self-interacting case (A = 0) and the self-interacting case (A = 1), without and 
with mass renormalization. The two horizontal bands indicate the values of 
exp(— Tim 2 /eE) for m 2 = 0.01 and m 2 = 0.27 (see the text for an explanation 
of these values). 

This problem is illustrated in the figure 13. The black curve shows the 
spectrum of produced particles (at t = 10, i.e. shortly after the external field 
has been switched on) without self-interactions (i.e. with a coupling constant 
A = 0). The (bare) mass in the Lagrangian, and hence in the classical equation 
of motion, is m = 0.1. This curve should be compared to the blue curve, where 
the same computation has been performed with a nonzero self-coupling A = 1, 
and the same value of the bare mass. We see that the particle yield has been 
considerably reduced. As we shall argue, this is a consequence of the fact that 
these two computations correspond to two different values of the renormalized 
mass. This is an unphysical effect that should be fixed. Indeed, one expects 
that the self-interactions among the scalar fields alter their long time evolution 
(and in particular play a crucial role in their thermalization), but should have 
little physical effects at short times. This issue is also visible in the two plots 
of the figure 14, where the computation at A = and the bare computation at 
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A = 1 lead to very different results, even at short times. 




Figure 14: Time evolution of the electrical current density (left) and of the 
electrical field (right). Curves are shown for the non self-interacting case (A = 0), 
and for the self-interacting case (A = 1), without and with mass renormalization. 

As explained in the appendix A, one can remove the quadratic divergence 
coming from the tadpoles by adding a mass countcrtcrm dm 2 in the equation 
of motion for the classical scalar field. This mass counterterm is space-time 
independent, and can thus be computed once for all at the initial time. We 
define it by 

Sm 2 = -A (ip* (0, x)(p(Q, x)) , (43) 

which is directly given by eq. (14) and is obviously independent of the external 
electrical field. With the value of the self-coupling A = 1 that we are using, we 
have A (ip*(0,x)(p(0,x)) — 0.26, which is why the calculation done without any 
mass renormalization gives a yield that is well reproduced by exp(— irm 2 /eE) 
with m 2 = 0.26 + 0.01 = 0.27 (see the figure 13). In the figure 13, we also 
show the p z spectrum obtained when this counterterm is added to the classical 
equation of motion. Now, we see that the particle yield is back at the level 
expected for a renormalized mass 14 m 2 R = 0.01. Similarly, the figure 14 shows 
that this mass renormalization cures the unphysical effect of the self-interactions 
at short times. 



6 Conclusions 

Motivated by recent works on thermalization in heavy ion collisions using clas- 
sical statistical field theory, in which one computes 1-loop quantum corrections 
by performing a Gaussian average over the initial condition of a purely classi- 
cal field, we have applied the same method to the calculation of the Schwinger 
mechanism of particle production in an external electrical field. 

14 The combination m 2 + Sm 2 that appears in the equation of motion should now be viewed 
as the bare mass. This bare mass combines with the tadpole that result from the Gaussian 
average, in such a way that the renormalized mass is back at the expected value. 
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We have first shown that, at leading order (i.e. at 1-loop), the spectrum 
of particles produced by the Schwinger mechanism can be expressed as a path 
integral over classical fields that have a Gaussian ensemble of initial conditions. 
The 2-point correlation function that characterizes this Gaussian distribution 
is uniquely determined by the propagator of small fluctuations on top of the 
external field. This representation of the spectrum is exact at 1-loop, and is 
a mere rewriting of the original quantum field theory result. Moreover, this 
formulation of the spectrum leads to a very efficient method for the numerical 
evaluation of the spectrum on a lattice. 

Then, by promoting the gauge potential to a dynamical variable, this formu- 
lation makes easy the inclusion of the back-reaction effects that are important 
for the physical consistency of the model. Indeed, the charged particles that 
are produced via the Schwinger mechanism tend to screen the applied external 
field, thereby reducing progressively their production rate. This back-reaction is 
essential for the conservation of total energy (particles + electromagnetic field). 

In the last section, we have studied the possibility of taking into account 
the self-interactions of the produce charged particles, which is an essential in- 
gredient for their eventual thermalization. In the classical statistical approach, 
they can be simply included by keeping the non-linear interaction term in the 
classical equation of motion for the field. However, since the Schwinger mecha- 
nism is very sensitive to the value of the mass of the particles, it is important 
to take proper steps in order to renormalize the mass: the naive (i.e. without 
any mass renormalization) inclusion of the self-interactions leads to an unphysi- 
cal -lattice spacing dependent- reduction of the charged particle yield, because 
the self-interactions produce large corrections to the mass. The main correc- 
tion to the mass is a quadratic ultraviolet divergence that comes from tadpole 
corrections. We have shown that it can be systematically subtracted in the clas- 
sical statistical framework by adding a mass counterterm in the classical field 
equation of motion. 

To close this paper, we would like to digress with a remark regarding the ap- 
plicability of the classical statistical approach to the calculation of observables. 
The example considered in this paper, as well as other quantities previously con- 
sidered in the literature, is an inclusive observable. This means that it measures 
the expectation value of a certain operator (here the particle number operator) 
in the final state of the system, without putting any constraints on this final 
state. This is the reason why these observables can be expressed in terms of 
fields that obey retarded boundary conditions, which in turn are amenable to a 
computation in terms of a Gaussian average over their initial conditions. Very 
little is known about exclusive observables, whose definition veto certain final 
states. The archetypal example of this kind of observable would be the probabil- 
ity of producing a specific number of charged particles, which obviously excludes 
any final state that does not contain the expected number of particles. In some 
examples, it has been shown that exclusive observables can be expressed at 
leading order in terms of classical fields that obey non-retarded boundary con- 
ditions (e.g. fields that are constrained both at t = — oo and at t = +00). At 
least on the surface, it seems very implausible that such an observable can be 
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obtained by a classical statistical simulation where one performs an average over 
the initial conditions of the classical field. 
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A Tadpoles in classical statistical field theory 

In this appendix, we show how to deal with the quadratic divergences that arise 
from the tadpoles when we keep the self-interactions in a classical statistical sim- 
ulation. In this appendix, we disregard the coupling of the scalars to the gauge 
fields, since our purpose is to discuss an issue related to the scalar self-coupling. 
Moreover, in order to keep the notations simple, we use continuum notations 
in this appendix. In practical applications to an actual lattice computation, all 
the integrals would become discrete sums. 

Let us illustrate the issue in the case of the simple calculation of the expec- 
tation value of the field operator itself, Before we go further, it is useful 
to recall the Green's formula that relates the classical field ip(x) to its initial 
value at x° = 0: 

= d 4 yG Jx 1 y)^(y) V 2 (y)+ f d 3 y G° R [x, y) «9>(0, y) , (44) 

L Jy°>0 Jy°=0 

where G° (x, y) is the free retarded Green's function of the Dalembertian opera- 
tor. By iterating the interaction term, one sees that the functional dependence 
of tp with respect to its initial condition can be represented as a sum of tree 
diagrams, of the form: 




In this equation, we have represented the terms that arise up to the second 
order in the coupling. The vertical dotted line symbolizes the initial time sur- 
face y° = 0, and the red circles the initial value of the field ip at this initial 
time. The average over the initial field is a Gaussian average; it can be done 
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diagrammatically by introducing the following objects 



<*> = * (!) = c 



(46) 



where the green dot represents the central value of the Gaussian ensemble 15 , and 
the green link represents its 2-point correlation function. When we apply these 
diagrammatic rules to the expectation value of the held operator, we obtain the 
following contributions 



In this diagrammatic representation, the green lines and dots represent the 
average over the initial condition, while the black lines are genuine (retarded) 
propagators coming from the subsequent time evolution. The first line is the 
sum of tree level contributions, the second line is the 1-loop contribution, etc.. 
The tree level contribution (first line) is nothing but the classical held whose 
initial condition at x° = is the central value of the Gaussian ensemble, ipo- We 
see that in the classical statistical approach, this classical solution is corrected 
by an infinite set of loop corrections. 

The tadpoles are the loops that have the strongest dependence on the ultra- 
violet cutoff. Before going further, it is worth clarifying an important point: it 
is not obvious a priori that the tadpole-like subgraphs that appear in the dia- 
grams of eq. (47) are identical to the usual tadpoles encountered in Feynman's 
perturbation theory. Let us demonstrate that they are in fact equal. By using 

15 Assuming a non-zero central value is a slight generalization of eq. (14), that allows to 
have a non-zero (4>(x)) for the purposes of the discussion in this appendix. 





(47) 
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the Green's formula (44), the tadpoles of eq. (47) can be expressed as 



d 3 yd 3 z G° R (x,y)d° y G° R (x,z)d° z (^(0, v)^(0,*)> 



y°,z°=0 



d 3 yd 3 z 



d^pd^q 



1 



(2tt) 8 {p 2 + ip°e)(q 2 +iq°e) 



y°,z°=0 



e -ip-(x-y) qO^ £ -iq-(x-z) gO _ 



d 3 k 
(2tt) 3 2/c 



(<p(0,y)<p(0,z)) 

(48) 



In the second line, we have replaced the retarded propagators by their Fourier 
representation. Then, a straightforward calculation, using eqs. (14), leads to 
the final expression - which is indeed the usual vacuum tadpole. Given the 
combinatorics for connecting a ip and a tp* in the product ip*<p 2 , each tadpole 
will arise with a prefactor A in the expansion of eq. (47) , which is also the right 
coupling and symmetry factor. 

It is possible to subtract all the tadpoles that arise during the time evolution 
by adding a mass counterterm Sm 2 in the classical equation of motion, that 
becomes 

(□ + m 2 + 5m 2 ) ip+ -<p*<p 2 = . (49) 

The effect of this counterterm in the equation of motion is to insert in its solution 
a mass counterterm in every place where a tadpole is allowed to appear. It 
therefore adds the following contributions to eq. (47) 



S (<p(x)) 




(50) 




where the red cross denotes the mass counterterm. The outcome is that the tad- 
poles are systematically cancelled by this procedure if one tunes the counterterm 
to precisely cancel the quadratic divergence of the tadpole, 



Sm 2 + A 



d 3 k 
(27r) 3 2fc 



< oo 



(51) 



Eq. (51) only specifies the divergent part of the mass counterterm; it also has 
a finite part that should be adjusted in order to have the desired renormalizcd 
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mass. Note also that the inclusion of this counterterm in the equation of motion 
only subtracts the quadratic divergences, possibly leaving a residual logarithmic 
dependence on the lattice spacing. 
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